clear, clc;

img = imread("haze_imgs\sam_3.bmp");
img = im2double(img);
[h,w,c]= size(img);
dark_channel = zeros(h,w);

%for i in range(1,10,1)(start,end,step)
%for i=start:step:end
for i=1:h
    for j=1:w
        % 求取暗通道
        dark_channel(i,j) = min(img(i,j,:));
    end
end
filter_d_c = ordfilt2(dark_channel,1,ones(3,3),'symmetric');
% filter_d_c = imerode(dark_channel,ones(3,3));
% t = 1 - I/A;
% 求A: 暗通道前0.1%最大值的值
top_num = max(1,round(h*w*0.001));
[sorted_dark, sorted_indices] = sort(filter_d_c(:),'descend');
top_indices = sorted_indices(1:top_num);
[rows,cols] = ind2sub([h,w],top_indices);
top_rgb = zeros(top_num,3);
for i = 1:top_num
    top_rgb(i,:) = img(rows(i),cols(i),:);
end
% 选择亮度最高的作为大气光值
[~,max_ind] = max(sum(top_rgb,2));
A = top_rgb(max_ind,:);
A = max(A);
% 求透射率, t = 1-w*(Id/A);
w0 = 0.95; t0 = 0.1;
t = 1-w0*(filter_d_c/A);
final_t = max(t,t0);
% 计算原图
out_pic = zeros(h,w,c);
for k=1:c
    for i=1:h
        for j=1:w
            out_pic(i,j,k) = (img(i,j,k)-A)/final_t(i,j) + A;
        end
    end
end
subplot(1,2,1);
imshow(img);
title("有雾图");
subplot(1,2,2);
imshow(out_pic);
title("传统DCP去雾");

